Realistic spatial and temporal earthquake distributions in a modified 

Olami-Feder-Christensen model 
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The Olami-Feder-Christensen model describes a limiting case of an elastic surface that slides on 
top of a substrate, and is one of the simplest models that display some features observed in actual 
seismicity patterns. However, temporal and spatial correlation of real earthquakes are not correctly 
described by this model in its original form. I propose and study a modified version of the model, 
that includes a mechanism of structural relaxation. With this modification, realistic features of 
seismicity emerge, that are not obtained with the original version, mainly: aftershocks that cluster 
spatially around the slip surface of the main shock and follow the Omori law, and averaged frictional 
properties similar to those observed in rock friction, in particular the velocity weakening effect. In 
addition, a Gutenberg-Richter (GR) law for the decaying of number of earthquakes with magnitude 
is obtained, with a decaying exponent within the range of experimentally observed values. Contrary 
to the original version of the model, a realistic value of the exponent appears without the necessity 
to fine tune any parameter. 

PACS numbers: 



I. INTRODUCTION 

The interest to describe the seismic phenomenon 
as originated in instabilities of dynamical systems has 
steadily increased in the last years. Although simplistic 
models cannot expect to reproduce the full phenomenol- 
ogy observed in real seismicity, it has become clear that a 
relatively fair understanding of some prominent features 
can be obtained using rather simple models. Histori- 
cally, the first model of this type is the one proposed 
by Burridge and Knopoff (BK) J- They considered a 
chain of elastically interacting rigid blocks (assumed to 
model portions of a tectonic plate) that are forced to 
slide onto an underlying surface. In order to obtain in- 
stabilities during motion that can be associated to earth- 
quakes, a crucial ingredient in this model is the use of 
a 'velocity weakening' friction force between blocks and 
substrate, i.e., a friction force that decreases as the rel- 
ative velocity increases. This was shown to generate in- 
stabilities that produce abrupt and potentially large re- 
arrangements of the blocks (the 'earthquakes') 0- The 
model is considered to be a paradigmatic case of self- 
organized criticality0, since the number of earth- 
quakes N within a fixed magnitude interval decays (al- 
beit in a limited magnitude range) exponentially with 
the magnitude Ad of the events, reproducing the empir- 
ical Gutenberg-Richter lawQ, namely N{M) ^ 10"''^^. 
For actual earthquakes the exponent 6 is usually found 
to be close to 1. 

A number of modifications and generalizations have 
been proposed to this model over the years. I will 
concentrate in the work of Olami, Feder, and Chris- 
tensen (0FC)3 who proposed a cellular automaton 
model based on the BK model, that has quite remark- 
able properties and is simple enough to be simulated very 
efficiently The OFC model considers a set of real val- 
ued variables Ui where i indicates the position in a two 
dimensional lattice. Ui is the force that the substrate ex- 



erts on the block at position i, and it represents the local 
stress between the sliding plates. The system is driven by 
uniformly increasing the values of Ui with time at a rate 
V , simulating the tectonic loading of the plates. Every 
time one of the variables Ui reaches a maximum pinning 
force (ordinarily set to an uniform, dimensionless value 
of 1), the local stress Ui is relaxed by setting it to zero 
(thus the local stress drop Au is equal to 1). The local 
stress drop produces a stress increase onto neighbor sites 
according to Uj — >■ Uj + a, where j indicates a neigh- 
bor site to i. The value of a can vary between and 
Uc = 1/ z, z being the number of neighbors in the lattice. 
The case a = ac is called the conservative case, whereas 
a < Uc are non-conscrvativc cases. A discharge event 
can produce the overpassing of the maximum local stress 
on one or more than one neighbors, and in this way a 
large cascade can be generated. This cascade is called 
an event, and is identified with an individual earthquake 
(note that the full cascade is assumed to occur at con- 
stant time, namely, earthquakes are instantaneous) . The 
size S of an event is defined as the sum of all discharges 
Au that compose the event, and the magnitude is de- 
fined as A/ = I logj^o 'S'l so to match (up to an additive 
constant) the usual definition used in geophysics Q. 

The OFC model is typically simulated using open 
boundary conditions. In the case of periodic boundary 
conditions the model exhibits a strong global synchro- 
nization originated in the spatial homogeneity. The OFC 
model displays an exponential decay of number of events 
as a function of magnitude compatible with a GR law. 
The 6-value is however not universal, but depends on 
the value of a. Realistic values of h are obtained for 
a ~ 0.2 (with z — A). A cut-off at large event sizes ex- 
ists that moves progressively to larger values for larger 
system sizes. 

After its introduction, the OFC model has been stud- 
ied in great detail, trying to extract from it the charac- 
teristics that are observed in actual seismicity patterns. 



2 



Although the finding of a GR decay law is a goal of this 
kind of model, the spatial and temporal clustering of 
earthquakes observed in real seismicity are certainly not 
reproduced by the OFC model (see below the discussion 
about aftershocks in the OFC model), as well as they 
were not reproduced neither by the model of BK. I refer 
in particular to the phenomenon of aftershocks, that has 
a partial quantitative description through the empirical 
Omori lawd This law states that the number of 

earthquakes in excess of its average value after a large 
event decays as ~ l/(t + c)P, where t is the time from the 
main shock, p is an empirical exponent, and c is a time 
constant in the range between minutes and hours. Usu- 
ally the value of p is found to be close to one, although 
other values and even other functional forms have also 



been proposed [11 



My contribution here is to modify the original OFC 
model in a way that allows for the existence of some 
kind of structural relaxation [l^. or aging in the sys- 
tem. This modification produces the appearance of cor- 
related events in the dynamical evolution, in particular 
aftershocks, generating earthquake sequences that qual- 
itatively and quantitatively resemble real ones. In addi- 
tion, the modified model will be shown to posses average 
friction properties that are compatible with those exper- 
imentally observed in rock friction studies [isj. In partic- 
ular, I obtain the effect known as 'velocity-weakening', 
namely, a reduction of the average friction force when 
the sliding velocity is increased, that is known to occur 
in rock frictioUjand plays a crucial effect in the triggering 
of earthquakes [g . Velocity weakening has been described 
phenomenologically in terms of the so called rate-and- 
state equations [l3| 1 but no detailed quantitative theory 
exists for it. 

In the next Section I introduce and justify the modifi- 
cations that are made onto the OFC model. Results arc 
presented in Section III. In Section IV, I show the de- 
pendence of some of the results on the kind of relaxation 
mechanism used. Finally, In Section V, I give some qual- 
itative interpretation of the appearance of aftershocks, 
summarize and conclude. 



II. THE MODIFICATIONS TO THE MODEL, 
AND THEIR JUSTIFICATION 



Two modifications will be implemented onto the origi- 
nal OFC model. They are the existence of random thresh- 
olds, and structural relaxation. I now present and justify 
them separately. 

1) Random thresholds: In the OFC model, the maxi- 
mum values that the variables Ui can withstand are set 
to a constant value of 1. Having in mind a realistic situa- 
tion of a heterogeneous fault, with the constitutive mate- 
rials having different properties at different positions, it 
becomes natural to consider a case in which the thresh- 
old values are not constant but have some spatial varia- 
tion. In concrete, the values of the local thresholds will 



be called uf^, and I draw them from a Gaussian distri- 
bution centered at uq, with standard deviation a. Each 
time Ui overpasses the local threshold u*'', Ui is updated 
to a new value. In concrete, I will use the update rule 
Ui —5- Ui — 1, i.e, I maintain (as in the original model) 
the prescription of a unitary local stress drop. Upon this 
drop of the local stress, the values of u on neighbor sites 
arc updated as before, namely Uj — > Uj+a, for j neighbor 
to i. 

Every time Ui is updated, a new value is assigned to the 
local threshold uf^ , taken from its original Gaussian dis- 
tribution. This prescription is justified on the same phys- 
ical arguments as before, since the sliding pieces can rea- 
sonably be thought to find different maximum strengths 
as sliding proceeds. I found that even a small value of a 
(about O.OSuq) is sufficient to qualitatively modify the be- 
havior of the OFC model (see results below). This means 
that the OFC model is not robust with respect to this 
perturbation, and since random variation of parameters 
is experimentally expected, the OFC model can prob- 
ably be considered as an interesting dynamical system, 
but not as a realistic model of the seismic phenomenon. 
In the language of renormalization group theory, we can 
say that fluctuations in the threshold values are "rele- 
vant" variables [Tsj . 

2) Structural relaxation: The existence of internal tem- 
poral effects in sliding systems is well established. Di- 
eterich and Kilgore[16j were the first to provide direct 
evidence that solid bodies in contact experience plas- 
tic relaxation that induces the increase of real contact 
area over time. In turn, this increase of contact area is 
intimately related to the velocity weakening effect and 
thus it affects directly the characteristics of the seismic 
phenomenon [sl. \l% . Thus the second modification I make 
on the OFC model is the inclusion of these relaxational 
processes [III . In order to justify the particular way in 
which relaxation will be introduced, I note that the plas- 
tic processes I am trying to model, always produce a 
reduction of the total energy E stored in the system. 
This energy will be dependent of the values of the Ui, 
namely there will be some functional form E[ui). The 
proposed relaxation mechanism causes a progressive re- 
duction of this energy through a standard first order re- 
laxation equation of the formpoj 



dui 
~dt 



R V 



,5E 



(1) 



where is the discrete Laplacian on the underlying 
square lattice, i.e., {V'^ f)i — fj ~^fi^ where j stands 
for the four neighbor sites to site i, and lattice parameter 
is taken as the unit of length (see also the discussion in 
Section IV). 

In order to have a concrete realization, we need to spec- 
ify the form of the function E{ui). Whereas in principle 
this can be done by a detailed derivation of the OFC 
model from an elastic model, at present I will take the 
view of proposing the simplest form for the relaxation 
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equation. This is obtained by taking E (x uf, and it 
gives 

^ = i?(vH + ^ (2) 

where the last term comes from the external driving, and 
a constant has been absorbed in the value of R. In the 
last section I will discuss the possibility of other analyti- 
cal forms of the relaxation mechanism and the effect on 
the results obtained. 

This mechanism[l2, produces (for V ^ 0) the pro- 
gressive uniformization of the local forces Ui on a time 
scale set by the relaxation parameter R. Thus, the rele- 
vant parameter of the dynamics of the system will be the 
ratio R/V, that measures the competing effect between 
relaxation and the global driving. 

Since I found that in the presence of non-uniform 
thresholds and for sufficiently large system sizes the re- 
sults become independent of the boundary conditions 
used, I chose to work always with periodic boundary con- 
ditions to reduce size effects as much as possible. Note 
also that the dynamics of the model is independent of 
the average value of the thresholds uq, since a change in 
this value produces only a rigid change of all Ui. I will 
formally take uo = 1, having in mind that a different 
value of Uq can be considered if we want, for instance, to 
maintain all Ui variables to be positive at all times, as its 
physical interpretation would suggest. In addition, very 
small events (those that span ten or less lattice points) 
are systematically cut off from the results, since they are 
spuriously dependent on the underlying numerical lat- 
tice. 



III. RESULTS 

Even in the case R = 0, there are qualitative dif- 
ferences between the results obtained with the present 
model (that uses random thresholds) and with the OFC 
model. In Fig. [T]we see results for this case {R — 0), 
for different values of a and a. The decaying exponent 
of the number of events with magnitude [Fig. [Ha)] is 
b ~ 0.37, independently of the precise values of a and 
a. There is a cut-off at large event sizes that moves to- 
wards infinity as a — ^ ac = 0.25. This is in contrast with 
the results in the OFC model, where the b exponent de- 
pends on a, and the cut-off depends on system size, and 
suggests that the physics of the model presented here is 
very different from that of the OFC model. Actually, the 
behavior I find for i? = is consistent with the case of 
an elastic interface driven on top of a disordered pinning 
potentialim. In particular the value of & ~ 0.37 com- 
pares very well with the avalanche size exponent for an 
elastic interface r ~ 1.25 (as defined for instance in [23l |. 
note that r = 1 + 2&/3). 

An inspection of the spatial and temporal sequences 
of epicenters (i.e., the triggering positions) of the events 
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FIG. 1: (Color online) (a) Magnitude-frequency distribution 
in the absence of relaxation (i? = 0) , and for different param- 
eters, as indicated. Continuous lines has a slope b = 0.37, 
for reference. The large size cut-off is mainly controlled by 
a, moving to infinity for a Oc = 0.25. There are no im- 
portant finite size effects in these results, as the system size 
{L X L = 200x200) is much larger than the largest event that 
is observed to occur for each value of a. (b)Magnitude vs 
time plot, and (c)projected position (along x-axis) vs time of 
epicenters of events with Af > 0.9, for the case for a = 0.2499, 
a = 1. Symbols size and color are magnitude dependent, ac- 
cording to the legend. No obvious sign of temporal or spatial 
correlation is observed. 

presented in Fig. [TJb)-(c), reveals no obvious correlations 
of any type. The conclusion from here is that the model 
with random thresholds and without relaxation {R = 0) 
is qualitatively different from the OFC model, but also 
far from being realistic in simulating seismicity. 

Qualitatively new results appear when R is different 
from zero. The magnitude frequency distribution for in- 
creasing values of R/V is shown in Fig. [2][a) for a fixed 
value of a = 0.246. As R/V increases, the b value in- 
creases [Fig. [2][a) inset]. Most remarkably, b seems to 
reach a well defined value (& ~ 1.0) when R/V is large, 
that is independent of a and <t (see Fig. [2ljb)) and that 
is comparable to actual values observed in earthquakes. 
Fig. [2{b) also shows that the large size cutoff of the GR 
plot is mainly governed by the value of a, moving to- 
wards infinity as a — ^ 0.25. The appearance of a realistic 
b value is encouraging since it is obtained without any 
tuning of parameters in the model (beyond the fact of 
R/V being sufficiently large). 
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FIG. 2: (Color online) (a) Magnitude-frequency distribution 
for increasing relaxation, and for a = 0.246, cr = 1, and 
L — 200. In the inset, the b values extracted from the main 
panel are plotted as a function of R/V, showing a convergence 
to a well defined value for large R/V. Panel (b) shows the 
independence of the b value on a and a, and the progressive 
increase of the large size cutoff as a — >■ ac = 0.25. 




FIG. 3: (Color online) (a) Magnitude-time plot and (b) pro- 
jected position (along x-axis) vs time of epicenters of events 
with M > 1.5, for the case R/V = 20, a = 0.246, a = 1. 
Temporal and spatial clustering is apparent. 



In addition to change the 6-exponcnt of the magnitude- 
frequency distribution, relaxation generates non-trivial 
correlations in the spatial and temporal event distribu- 
tion. We will see now that this clustering has features 
that are known to correspond to real earthquakes. First 
of all, it is qualitatively seen in Fig. [3] that events ac- 
cumulate following large ones, and that the epicenters 
of the clustered events occur close to the epicenter of 
the large shock. This reproduces well known features 



of real aftershocks. In order to provide a more quan- 
titative characterization of aftershocks in the model, in 
Fig. HI a), I show the result of calculate histograms of 
events occurrence around main shocks. For this analysis, 
a main shock is operationally defined as any event having 
M > 3.0 (this corresponds to events producing a rupture 
region of linear size about 60 lattice sites). Time is set to 
zero at the main shock. Different curves are presented, 
that correspond to events occurring within a given spa- 
tial distance d from the main shock epicenter. Curves are 
normalized in such a way that N{t — )■ oo) = 1. The over 
abundance of events following large ones is clear. There 
is also some over abundance of events preceding large 
ones (foreshocks), but in a much lesser extent than the 
case of aftershocks. In order to compare with an Omori 
expression, in Fig. IDJb) I plot the evolution of the activ- 
ity after the main shock, in logarithmic scale, and using 
different lower cutoff values to define aftershocks. Curves 
are compatible with an Omori law, but we see that the 
activity soon becomes masked by the background activ- 
ity, rendering a quantitative determination of parameters 
in the Omori expression very unconstrained, (below we 
will see a trick to avoid this problem) . In addition, we see 
that the convergence towards the background activity is 
not monotonous, instead a time window of lower-than- 
average activity (at times around 0.05/i?) is observed. 

In order to analyze quantitatively the spatial cluster- 
ing of aftershocks, in Fig. |4l^c) I present histograms of 
activity in the system for different time windows after 
the main shock as a function of the distance to the main 
shock epicenter. It is observed that in very short times 
after the main shock, aftershocks occur rather uniformly 
in a region of size comparable with the rupture region of 
the main shock, and fewer events are observed at larger 
distances. When we consider later aftershocks, the spa- 
tial distribution clearly drifts away of the main shock 
epicenter. This is in nice agreement with the observed 
behavior of actual seismicity patterns. The lower-than- 
average activity region indicated in (b) is seen here [latest 
curve in (c)] to be due to the eventual appearance of a 
region of depleted activity (with respect to its spatial 
average) close to the main shock epicenter. 

A complementary, more visual example of aftershocks 
spatial distribution is presented in Fig. [S] There I show 
the actual region that was broken by a particular main 
shock, and the activity following this event. Squares, 
circles, and triangles correspond to three time windows of 
progressively later events, as indicated in the legend {t = 
at the main shock). Size of the symbol increases with 
magnitude. In this example we see again that aftershock 
activity starts mainly within the region broken by the 
main shock, and then progressively drifts away. 

In the previous analysis, the aftershock statistics is lim- 
ited by the fact that tectonic loading continues to trigger 
events in the system that rapidly mask the aftershock 
tail of previous main shocks. However a numerical trick 
can be implement to overcome this problem in order to 
study aftershocks in more detail. Along a simulation. 
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FIG. 4: (Color online) (a) Aftershocks and foreshocks in 
the model with parameters R/V = 20, a = 0.244, a — 1, 
system size 200 x 200. I show the accumulated histogram 
N{t) of events occurring at time t with respect to the large 
event (I consider about 200 large events, with magnitude 
M > 3.0, with relative spatial distance of epicenters smaller 
than the cut off value d . The histogram is normalized to 
the corresponding totally random situation [this means that 
N{t — >■ oo) = 1]. The accumulation of events following the 
large ones (aftershocks) is clearly visible. Foreshocks are ob- 
served but in a much lesser extent, (b) Activity in the whole 
system following main shocks, in logarithmic scale. For com- 
parison, an Omori expression a/{t + cY (with a — 0.035, 
c = 0.001, p — 0.8) is also plotted, (c) Activity after main 
events, as a function of the distance to main shock epicenter d, 
for different time windows after the main event, as indicated. 
The typical radius of regions broken by the main shock is indi- 
cated by the vertical dotted line. A drift away from the main 
shock epicenter of later time aftershocks is clearly visible. 
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FIG. 5: (Color online) In gray, the broken region of a partic- 
ular main shock. Symbols are the epicenters of events follow- 
ing the main shock. Symbol size depends on magnitude, and 
symbol type indicates different time windows after the main 
shock. The figure displays a system of size 200 x 200, and 
parameters are a = 0.244, R/V — 20, and a — 1. 




FIG. 6: (Color online) Aftershock decay rate in simulations 
in which tectonic loading is stopped {V = 0) after a main 
shock has occurred (parameters as in Fig. |4|. The asymp- 
totic form of the decay follows nicely an Omori decay (for 
reference, the continuous lines has a slope of 1.1). Different 
curves correspond to different lower cutoff values Mq used to 
count aftershocks. 



once a main shock is detected I stop tectonic loading 
(i.e. setting V = 0) and look for the aftershock occur- 
rence. The process is repeated many times to accumulate 
good statistics. The result is shown in Fig. [HI Now we 
can follow the aftershock decay rate for few orders of 
magnitude in time. It is seen that a very good power law 
(Omori like) decay is obtained with a p exponent around 
1.1. I notice however that at short times after the main 
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FIG. 7: (Color online) (a) Average friction force _F as a 
function of velocity in a system of size 200 x 200 sites and 
other parameters as indicated. Velocity weakening follow- 
ing an approximately logarithmic dependence on velocity is 
clearly observed, (b) Results of three individual realizations 
for the stress decay in a system of size 200 x 200 (R/V = 20, 
a — 0.244, (7 = 1) after the abrupt stopping of loading at 
some arbitrary time (set to zero). 



shock some over abundance of events with respect to the 
asymptotic rate is observed. This can be interpreted as 
a larger p- value if the aftershock sequence is observed in 
a limited time interval. 

It is necessary to mention here that aftershocks of a 
very peculiar type have been found in the original OFC 
modeljlll- In my view, these aftershocks reflect once 
more the kind of synchronization that the OFC model 
is prone to, and have nothing to do with aftershocks ob- 
served in real seismicity. In particular, they completely 
disappear once a randomness in the thresholds of about 
5% is included [2^, indicating that they are not a robust 
property. Also, their existence depends exclusively on 
external loading, i.e., if at some moment the external 
driving vanishes, all seismic activity ceases immediately. 
Namely, the plot equivalent to that in Fig. |6]for the OFC 
model would be completely void. 

In addition to generate a realistic b value, and pro- 
duce clustering effects, structural relaxation generates 



global frictional properties that are comparable to what 
is observed in laboratory experiments of friction between 
solids. I refer in particular to the so-called velocity weak- 
ening properties of the friction phenomenon^, and 
in general to the phenomenology given by rate-and-state 
equations [T^. widely used in seismological analysis. Ve- 
locity weakening means that the average friction force F 
between the sliding bodies decreases as a function of rel- 
ative sliding velocity. In the BK model, this behavior has 
to be introduced by hand in the form of a tailored fric- 
tion law between blocks and substrate. The OFC model, 
on the other hand, can be considered to generate a fric- 
tion force (that is obtained in this case as the spatial and 
temporal average of the local friction forces Ui) that is in- 
dependent of sliding velocity, since sliding velocity plays 
no role in the dynamics of the system, as earthquakes are 
assumed to be instantaneous. But in the present model, 
the interplay between structural relaxation and the ex- 
ternal driving velocity generates a friction force F that 
depends on velocity. This dependence is of the velocity- 
weakening type, as can be seen in the results in Fig. ^a). 
The decrease of friction force with velocity is mainly log- 
arithmic in about three orders of magnitude of velocity 
variation, quite comparable to the experimental results 
m Ref. hj. The increase of friction force as velocity is 
reduced can be qualitatively understood if we consider 
that at lower velocities, the system has more time to 
reach more stable configurations (with lower energies). 
Thus, larger forces have to be applied in order to, even- 
tually, take out the system from these configurations, and 
this means a larger friction force. 

Another effect that the present model reproduces is the 
stress relaxation that is observed after loading is stopped 
in laboratory friction experiments [l^ . We have already 
seen that in the present model activity decays slowly after 
tectonic loading is stopped, and this produces a stress 
relaxation that follows an almost logarithmic trend. In 
fact, in Fig. El^b) we see examples of the stress decay 
after stopping loading, where the effect of individual large 
aftershocks is seen as abrupt stress drops. 



IV. DEPENDENCES ON THE RELAXATION 
MECHANISM 

The use of a conserving form of Eq. ((T]), i.e., the 
inclusion of the Laplacian operator, instead of a non- 
conserving form of the kind du/dt = —\5E/5u is difficult 
to justify on first principles. One possible a posteriori 
justification is that a non-conserving dynamics produces, 
in the absence of tectonic loading, the evolution of the 
system towards a state with Ui = everywhere, i.e., a 
stress free state. In particular, the almost logarithmic 
stress decay obtained in the last Section would not oc- 
cur, instead we would observe an exponential decay to- 
wards zero stress. This is not a realistic situation for the 
present problem, although it could possibly be appro- 
priate to model a viscoelastic response. The necessary 
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FIG. 8: (Color online) Magnitude frequency distribution ob- 
tained with the three different relaxation mechanisms pro- 
posed in text, for R/V — 20, a = 0.246, and cr = 1 in all cases 
(curves are not normalized, and were vertically displaced for 
better comparison). The value of b for the three mechanisms 
is within a range of 0.9 < b < 1.1, close to experimentally 
observed values. The two limiting slopes are plotted as con- 
tinuous lines, for reference. 



FIG. 9: (Color online) Time decay of aftershocks in a 200 x 
200 system in the case in which loading is stopped after a 
main shock is detected. The three curves correspond to the 
three different realizations of the relaxation mechanism. The 
curve corresponding to Eq. (3) is plotted in a different scale 
in the inset, to highlight an exponential decay of aftershocks 
in this case. 



requirement to have a system that is able to maintain a 
constant stress under static conditions is that relaxation 
does not modify the spatial average u of u,;. Eq. ([2]) 
certainly satisfies this requirement, but other forms are 
possible. I compare in this Section some of the results 
obtained using Eq. ([2]) with two other possibilities for 
the relaxation mechanism, namely 



(3) 
(4) 



(I use the same symbol R for the relaxation parameter in 
all cases, although its numerical value may differ). Both 
Eqs. ([3]) and ([4]) do not modify the spatial average u of 
the Ui variables. Eq. ([3]) is a sort of "mean field" imple- 
mentation of the relaxation mechanism. It is not a realis- 
tic possibility in a physical system with local interactions, 
but is an interesting case of study to compare with. Eq. 
^ incorporates a double Laplacian to drive the temporal 
variations of w^. I will consider this to be simply another 
formal possibility, although it can be given a more phys- 
ical justification by deriving the present model as a limit 
of an elastic spr ing block-model in the presence of relax- 



ation (see Ref . [21| ) . Now I will present a comparison of 
the results obtained using Eqs. ([2]), ([3]), and (|4|). 

First of all, the three implementations produce a mod- 
ification of the b exponent of the GR decay, that tends to 
a conserved value when R/V is large enough, in the three 



cases. The conserved value, is within the range 0.9-1.1 
in the three cases. This can be seen in Fig. HI where 
I present the magnitude frequency distribution for the 
three mechanisms, for a rather large value of R/V (so 
the b value has already reached its asymptotic value). 

The other comparison refers to the asymptotic decay of 
aftershocks. In Fig. IHlwe observe the aftershock activity 
(using the prescription of stopping tectonic loading) after 
main shocks. There is a clear difference in aftershock 
behavior for the mean field case of Eq. ([3]) . In this case, 
the aftershocks activity decays exponentially with time. 
On the other hand, the difference between the behaviors 
using Eq. ^ and (0]) is more subtle. I already indicated 
that using Eq. [5] an over-abundance of aftershocks at 
small times is observed, before an Omori decay (with 
p ~ 1.1) is observed. For the bi-Laplacian relaxation the 
aftershock excess is much less pronounced, and an almost 
perfect Omori law with p ~ 1.25 is observed in this case. 

The results in this Section indicate that some quali- 
tative features of the model are robust upon a change 
of the relaxation mechanism, although differences in the 
details can be expected. 



V. FURTHER ANALYSIS AND CONCLUSIONS 

Finally, in order to clarify the origin of aftershocks in 
the model, we can make an analysis of the limiting case in 
which V/R ^ 0. This case can be realized in the follow- 
ing way. We set temporarily V to zero, and allow only 
the evolution given by the relaxation term in Eq. 
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This evolution is continued until we can guarantee that 
no other event will be triggered by the relaxation alone. 
At this point the values of Ui can be set everywhere equal 
to their mean value, and this flat interface can be driven 
by the external velocity V until a new instability occurs. 
In this limit, a precise distinction between main shocks 
and aftershocks can be given: aftershocks are events that 
are triggered by the term proportional to R in Eq. 
while main shocks are triggered by V (note that in this 
case, main shocks are not defined in terms of their inten- 
sity, in fact, it can occur that a main shock produces an 
aftershock of larger magnitude that the starting event). 
The mechanism of aftershock triggering is illustrated in 
this limit in Fig. [101 In order for the relaxation to be 
able to trigger events by itself, the thresholds cannot be 
uniform, since in that case, starting with Ui < 1 after 
a given event, evolution through Eq. ([2]) with V ^ 
cannot produce any > 1. However, if thresholds have 
some randomness, the evolution according to Eq. ([2]) 
with y = can produce Ui > u^^ at some position (par- 
ticularly at those with the smallest thresholds), and an 
aftershock is triggered. This highlights the crucial role 
played by a non-uniform distribution of thresholds in the 
appearance of aftershocks. It is thus not surprising that 
aftershocks are observed only if the distribution of thresh- 
olds has a dispersion a larger than some minimum value 
that turns out to be about 0.25. 

Summarizing, I have presented a model that is based 
on the one proposed by Olami, Feder, and ChristensenQ 
to study the dynamical appearance of slip events between 
tectonic plates. Modifications consider the existence of 
random thresholds, and the possibility of relaxation in 
the system. Relaxation acts trying to strengthen the 
contact between the sliding surfaces if they remain at 
rest relative to each other. When the surfaces slip, the 
contacts refreshen and there is a competence between 
the relaxation mechanism and the external driving of the 
system. With this kind of modification I have been able 
to generate earthquake sequences that contain many of 
the features observed in real scismicity. In particular, I 
have presented results of temporal clustering of events 
following a main shock according to the Omori law, and 
spatial clustering of these events around the epicenter of 
the main shock. Also, an appropriate decay of number 
of events as a function of magnitude according to the 
Gutenberg-Richter law, with a realistic 6-exponent has 
been obtained, and the value of h was shown to be in- 



dependent of other parameters of the model, relaxation 
assumed to be sufficiently large. Although the present 
model does not introduce velocity weakening directly, 
this effect appears as a consequence of structural relax- 
ation. In this way, the model gives also as a physical 
basis for the rate-and-state equations used to describe 
frictional properties of solids. 
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FIG. 10: (Color online) Pictorial description of main shocks 
and aftershocks in a one dimensional version of the model, in 
the limit V/R ^ Q. Vertical position of thick (thin) segments 
are the values of Ui {uf^). Rate of change of u's is indicated 
by the arrows. In (a), a relaxed (thus uniform) distribution 
of It's increases at a rate V . At the point indicated by A, the 
first threshold will be overpassed and a main shock will be 
triggered. After the main event, the system reaches a situa- 
tion depicted in (b) . Relaxation tends to make the values of u 
uniform across the sample, and in this process an aftershock 
(in this case at point B) can be triggered. 
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